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1. 


INTRODUCTION 


Stringent  new  requirements  to  reduce  gravity- induced 
uncertainties  in  inertial  systems,  to  control  Geodetic  and 
Geophysical  (G&G)  factors  during  the  test  and  evaluation  of 
inertial  instruments,  and  to  achieve  mobility  for  new  classes 
of  ballistic  missile  systems  motivate  more  accurate  means  to 
measure  the  earth's  gravity  field.  These  requirements  demand 
increasingly  fine  local  gravity  data,  more  accuracy  at  high 
spatial  frequencies,  and,  in  many  cases,  coverage  of  large 
areas . 


In  this  report,  high- resolution  digital  terrain  ele¬ 
vation  data  (DTED)  are  used  to  compute  the  effect  of  terrain 
variations  on  the  very  high-frequency  vertical  gravity  field 
in  a  local  area  (wavelengths  shorter  than  10  km) .  The  terrain 
effect  is  computed  with  a  specified  root-mean -square  (rms) 
accuracy  as  determined  by  an  error  covariance  analysis  that 
accounts  for  the  finite  extent  of  the  DTED. 

The  terrain  effect  is  only  a  portion  of  the  total 
gravity  disturbance,  because  regional,  continental,  and  global¬ 
sized  density  variations  deep  below  the  earth's  crust  also 
contribute  significantly.  Nevertheless,  at  very  short  wave¬ 
lengths  the  terrain  effect  is  expected  to  dominate  the  gravity 
field  in  most  areas  (Ref.  1).  Therefore,  the  high-frequency 
terrain  effect  is  used  as  a  reasonable  estimate  of  the  gravity 
disturbance  at  wavelengths  shorter  than  10  km. 

This  report  describes  a  new  analytical  technique  for 
using  DTED  to  compute  the  high-frequency  terrain  effect  on 


gravity.  The  methodology  is  applicable  for  estimating  all 
three  components  of  the  high-frequency  gravity  vector  from 
DTED,  for  determining  the  lowest  spatial  frequency  appropriate 
for  the  estimation,  for  sizing  the  zones  of  DTED  coverage  needed 
for  a  specified  accuracy,  and  for  applying  modern  vector/parallel 
numerical  processing  techniques  on  supercomputers.  The  approach 
described  here  is  novel  in  two  ways:  it  uses  a  direct  quadrature 
approach  to  compute  the  terrain  effect,  and  it  includes  systematic 
control  of  accuracy. 

While,  in  principle,  the  accuracy  of  DTED  could  be 
approximately  1  to  3  m  rms ,  in  practice  much  of  the  data  now  in 
place  are  of  unknown  accuracy.  Thus,  the  results  presented 
herein,  as  well  as  widespread  use  of  the  methodology  described 
in  the  sections  which  follow,  appeal  to  the  future.  This  future 
will  be  realized  when  DTED  accuracy  can  confidently  be  known 
to  match  its  current  precision. 


2. 


BACKGROUND 


The  problem  addressed  in  this  report  is  to  estimate 
the  high-frequency  vertical  gravity  field  of  the  earth  at 
points  on  a  horizontal  grid.  The  vertical  component  of  gravi¬ 
ty  is  to  be  estimated  at  points  that  are  uniformly  spaced  at 
three  arc  second  intervals  over  a  rectangular  area  containing 
several  hundred  thousand  such  points.  The  longest  wavelengths 
of  interest  are  about  10  km.  Because  these  wavelengths  are  much 
smaller  than  the  radius  of  the  earth,  the  problem  is  formulated 
using  the  flat-earth  approximation  (i.e.,  the  reference  ellipsoid 
is  approximated  locally  by  a  horizontal  plane). 


2 . 1  PROBLEM  STATEMENT 

The  technical  approach  taken  in  this  work  is  to  use  the 

local  terrain  effect  on  gravity,  suitably  filtered  to  attenuate 

low-frequency  components,  as  the  desired  gravity  estimate.  This 

local  terrain  effect  is  defined  as  the  gravity  due  to  the  earth's 

mass  above  a  reference  plane  in  a  local  geographical  area  that 

is  large  enough  to  provide  a  specified  rms  accuracy  at  high 

spatial  frequencies.  To  compute  the  terrain  effect,  Newton's 

law  of  gravitation  is  used  to  express  the  vertical  gravity  as  a 

volume  integral  over  the  local  mass.  The  mass  density  of  the 

3  3 

earth  is  modeled  as  constant  at  2.67  x  10  kg/m  ,  which  is  the 
mean  crustal  density,  and  the  surface  of  the  earth  is  represented 
by  DTED  on  a  uniform  grid.  The  DTED  are  expressed  as  heights 
above  sea  level  to  the  nearest  meter. 


The  terrain-effect  integral  is  formulated  as  illus¬ 
trated  in  Fig.  2.1-1,  which  depicts  a  vertical  cross-section  of 
the  earth's  crust.  A  Cartesian  coordinate  system  is  used,  with 
variables  x  and  y  measuring  position  in  the  horizontal  plane  and 
variable  z  measuring  height.  The  observation  point  P  has  coordi¬ 
nates  (Xp,yp,Zp);  this  is  a  generic  point  at  which  the  vertical 
gravity  is  to  be  computed.  The  reference  surface  is  at  height 
Zq,  and  the  earth's  surface  at  position  (x,y)  has  height  h(x,y). 
The  terrain  effect  at  point  P  is  the  integral  of  the  gravity 
contributions  from  all  infinitesimal  mass  elements  such  as  the 
one  at  point  (x,y,z)  in  Fig.  2.1-1.  The  gravity  potential  at 
point  P  due  to  all  mass  elements  is  denoted  «j>(P)  and  is  given 
by  Newton's  law: 


4>  <P) 


oo 

■  #vvvv-’-- 


(2.1-1) 


-11  3  -1  -2 

where  G  is  the  gravitational  constant  (6.672x10  m  • kg  *s  ), 

3 

p  is  the  mass  density,  which  is  modeled  as  a  constant  (2.67x10 
_  3 

kg*m  ),  and  is  distance  from  the  observation  point.  The 
vertical  component  of  the  terrain  effect  at  point  P  is  denoted 
g(P)  and  is  measured  positive  in  the  downward  (negative  z) 
direction.  From  potential  theory,  it  follows  that 


Therefore , 


g(P) 


94>(P) 

dZ 

P 


g(P) 


h(x,y) 


z0 


(2.1-2) 


dz 


(2.1-3) 


G-02023 


Figure  2il-l  Vertical  Cross-Section  of  Earth's  Crust 

in  Relation  to  Observation  Point  P 


The  integration  with  respect  to  z  in  Eq.  2.1-3  can  be 
expressed  in  closed  form.  For  mathematical  convenience,  this 
integration  is  performed  in  two  steps: 


The  decomposition  of  Eq.  2.1-4  leads  to  the  following  expression 
for  the  vertical  terrain  effect: 


(2.1-5) 


g(p)  =  27iGp|zp-zQ| 


*  "IP 

-00 


where  q  and  r  are  defined  by  Fig.  2.1-1  and  by  Eqs .  2.1-6,  2.1-7 
and  2.1-8: 


q(x  ,y  ,x,y) 

d(Zp,x,y) 


P-X)2  +  (yp-y) 
-  h(x,y) 


(2.1-6) 

(2.1-7) 

(2.1-8) 


In  Eq.  2.1-5,  the  terrain  effect  is  expressed  as  the  sum  of  two 
terms ,  which  are  termed  the  slab  effect  and  the  terrain  correc¬ 
tion  .  The  slab  effect  is  the  first  term  and,  analogously  to  a 
Bouguer  plate,  represents  the  vertical  component  of  gravity  due 
to  an  infinite  horizontal  slab  of  mass  filling  the  space  between 
the  lower  reference  surface  and  the  observation  point.  The  ter¬ 
rain  correction  is  th^  two-dimensional  integral  over  the  hori¬ 
zontal  plane,  which  accounts  for  terrain  deviations  from  the 
top  of  the  slab. 

In  this  analysis,  all  observation  points  have  the  same 
heights  so  that  zp  is  fixed,  and  the  slab  effect  is  constant. 
The  major  challenge  in  computing  the  terrain  effect  is  the 
numerical  evaluation  of  the  terrain-correction  integral  in 
Eq.  2.1-5,  i.e. 
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I(P)  =  Gp 


■j^"  (r  '  q)  **& 


(2.1-9) 


Integral  1(P)  must  be  evaluated  for  each  grid  point  P  at  which 
gravity  is  being  computed. 


2 . 2  PREVIOUS  APPROACHES 

This  section  reviews  previous  approaches  to  the  calcu¬ 
lation  of  terrain  corrections  on  a  dense  grid  using  DTED:  the 
prism  approach  and  two  variations  on  the  more  recent  Fast 
Fourier  Transform  (FFT)  approach. 


2.2.1  PRISM  Approach 

A  frequently  used  approach  to  evaluating  the  terrain 
effect  based  on  DTED  is  to  approximate  the  terrain  as  a  col¬ 
lection  of  prisms  of  constant  mass  density  and  then  apply  an 
analytic  expression  for  the  gravity  disturbance  due  to  each 
prism  (e.g.  Ref.  1).  Flat-topped  prisms  are  often  used  in 
which  h(x,y)  is  assumed  constant  on  each  of  a  set  of  rectangles 
in  the  x-y  plane.  The  formula  for  each  prism  requires  the  evalu¬ 
ation  of  eight  logarithms  and  eight  arctangents,  along  with 
numerous  subtractions.  This  approach  imposes  several  drawbacks: 
the  high  computational  cost  of  evaluating  the  transcendental 
functions,  ill-conditioning  due  to  cancellation  errors,  and 
the  introduction  of  spurious  high  frequency  gravity  signal  due 
to  the  abrupt  edges  of  the  prisms. 


A  recently  developed  alternative  is  the  FFT  approach 
(Refs.  2,  3,  and  4),  in  which  the  terrain-correction  integral 
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is  approximated  as  a  sum  of  linear  convolutions  that  are 
implemented  as  multiplications  in  the  spatial  frequency  domain. 

The  key  approximation  in  evaluating  Eq.  2.1-9  is  to  replace 

1  _  1 
r  ’  q 

by 

1  d2 

-2  ? 


Note  that  this  approximation  is  quite  accurate  for  q>>d.  Further 
discussion  and  application,  accounting  for  this  restriction, 
is  presented  in  Section  3.2.  The  result  is 


(2.2-1) 


2.2.2  FFT  Approach 


For  the  case  of  the  classical  terrain  correction,  in 
which  the  computation  point  P  is  on  the  surface  of  the  terrain, 
the  FFT  method  is  developed  with  the  simplification,  z  =  hp , 
where  hp  =  h(xp,yp)  and 


d  =  hp  -  h(x,y) 


(2.2-2) 


The  first  term  on  the  right  side  of  Eq.  2.1-5  then  reduces  to 
the  classical  Bouguer  reduction  (e.g.,  Ref.  5).  Note  that  for 
zp  /  hp  the  integral  in  Eq.  2.2-1  fails  to  converge.  (A  simple 
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example  of  this  divergence  is  the  trivial  case  in  which  zp  /  0 
and  h(x ,y )  =  0 . ) 


IT. 


The  next  step  in  the  FFT  method  is  to  define 


so  that 


f(x,y)  =  (x2  +  y2) 


-3/2 


~3  E  f<vx>  yp-y> 

q 


(2.2-3) 


(2.2-4) 


Equation  2.2-1  may  then  be  written  as 


I(P)  3 


00 

>// 


(h(x ,y ) -hp)  f(xp-x,yp-y)  dxdy 


(2.2-5) 


Expanding  the  squared  quantity  leads  to  a  sum  of  linear  convo- 

2 

lutions  of  f  with  h  and  h  ,  plus  a  constant  term.  The  convo¬ 
lutions  may  be  evaluated  (approximately)  by  multiplication  in 
the  discrete  spatial  frequency  domain  with  the  use  of  a  two- 
dimensional  FFT.  Clearly,  one  advantage  of  such  a  method  is 
that,  upon  return  from  the  frequency  domain,  the  values  of  the 
terrain  effect  are  available  on  an  entire  grid  of  points. 

One  troublesome  problem,  however,  is  that  the  integral 
of  f  over  the  x-y  plane  is  divergent,  which  causes  the  Fourier 
transform  of  f  to  become  undefined.  Furthermore,  expanding 
Eq.  2.2-5  results  in  a  difference  of  infinite  quantities.  None¬ 
theless,  according  to  Refs.  3  and  4,  replacing  f(0,0)  by  0  in 
the  discrete  computations  can  lead  to  reasonable  numerical 
results.  In  this  instance,  the  argument  is  made  that  the 
integral  of  f  over  the  plane  can  be  replaced  by  the  discrete 
Fourier  transform  of  the  discretized,  modified  f,  evaluated  at 
zero  frequency.  Unfortunately  the  process  is  ad  hoc  and  uni¬ 
versal  applicability  cannot  be  assured. 
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A  variation,  proposed  in  Refs.  2  and  4,  is  to  replace 


f  by 


“3/2 

g(x,y)  =  (xz  +  yz  +  bz) 


(2.2-6) 


where  b  is  a  nonzero  constant  whose  dimension  is  length.  It 
can  be  shown  that 


g(x,y)  dxdy  =  2n/b 


(2.2-7) 


and  that  the  Fourier  transform  of  g  exists  and  has  a  simple 
analytic  form. 

2.2.3  Modified  FFT  Approach 

References  2  and  4  advocate  choosing  as  small  a  value 
of  the  parameter  b  as  possible.  However,  Ref.  6  shows  that 
substitution  of 


2  2 
(q2  +  bZ) 


Rp  =  (q2  +  (zp-z)2) 


(2.2-8) 


in  the  rightmost  integral  of  Eq.  2.1-4  leads  directly  to  the 
approximation 


I(P)  = 


00 

'  1  Gp  j| 


(h(x,y)-Zp)  g(xp-x,  y  -y)  dxdy 


(2.2-9) 


2 

instead  of  Eq.  2.2-5.  Equation  2.2-9  shows  that  b  should  be  a 

2 

kind  of  average  of  (zp-z)  and  should  not  be  arbitrarily  chosen 

to  be  as  small  as  possible.  Note  also  that,  in  the  formulation 

of  Eq.  2.2-9,  the  FFT  approach  is  no  longer  limited  to  computation 

points  on  the  terrain  surface;  zp  in  Eq.  2.2-9  need  not  equal  hp , 

and  I(P)  can  be  readily  evaluated  in  terms  of  the  convolution 

2 

of  g  with  h  and  h  . 

However,  even  with  the  improvements  afforded  by 
Eqs.  2.2-6  and  2.2-9,  an  approximation  is  involved  in  the  FFT 
method,  which  (like  the  prism  representation)  may  incur  high- 
frequency  errors  and  edge-effect  problems.  Therefore,  instead 
of  applying  one  of  the  conventional  methods,  the  preferred 
technique  pioneered  herein  is  a  completely  alternative  approach 
to  the  computation  of  terrain  corrections.  It  is  simple  in 
principle  yet  appears  to  have  been  overlooked  previously  because 
of  the  unavailability  of  modern  parallel  processing  resources: 
direct  quadrature  of  Eq.  2.1-9. 


This  chapter  discusses  several  issues  related  to 
the  numerical  integration  in  Eq.  2.1-9: 


Avoiding  cancellation  errors  in  the  eval¬ 
uation  of  the  integrand  using  finite  pre¬ 
cision  arithmetic. 

Accounting  for  the  singularity  at  the 
origin  in  Eq.  2.1-9. 

Determining  the  extent  of  the  integration 
domain  needed  to  achieve  the  required 
accuracy  for  short-wavelength  gravity. 


The  next  three  sections  discuss  each  of  these  issues. 


3.1  AVOIDING  CANCELLATION  ERRORS 

Equation  2.1-9  is  unsuited  for  numerical  computations, 
because  the  integrand  is  the  difference  of  two  terms  nearly 
equal  to  each  other.  When  finite  precision  floating  point 
arithmetic  is  used,  the  mantissas  of  both  terms  almost  coincide, 
and  the  difference  between  the  two  numbers  is  reflected  only  in 
their  least  significant  bits.  The  relative  accuracy  with  which 
the  difference  can  be  computed  is  equal  to  that  fraction  of  the 
mantissa  where  the  two  terms  differ. 

The  errors  that  arise  from  computing  the  difference 
between  two  nearly  equal  numbers  are  usually  termed  cancella¬ 
tion  errors.  To  avoid  them,  Eq.  2.1-9  can  be  written  as 


riw.  yrmTBwviWb 


KP)  = 


00 

'// 


ff  ^ 


dxdy 


(3.1-1) 


which  reduces  to 


I(P)  =  -Gp 


qr(r+q) 


dxdy 


(3.1-2) 


The  integrand  in  Eq.  3.1-2,  although  more  complicated  than  the 
integrand  of  Eq.  2.1-9,  is  now  much  better  suited  for  numerical 
computations . 


3.2  ACCOUNTING  FOR  THE  SINGULARITY  AT  THE  COMPUTATION  POINT 


The  integrand  in  Eq.  3.1-2  is  singular  at  the  computa¬ 
tion  point  P,  where  (x,y)  =  (Xp.y^).  By  changing  to  polar  coor¬ 
dinates  centered  on  the  computation  point,  it  can  be  seen  that 
the  singularity  at  the  computation  point  is  integrable.  Thus, 
from  a  purely  theoretical  point  of  view,  the  singularity  pre¬ 
sents  no  problem.  However,  the  limiting  process  needed  to 
integrate  the  singularity  requires  the  knowledge  of  h(x,y) 
over  a  continuum  in  the  neighborhood  of  (x^.y^).  Since,  in 
practice,  the  terrain  elevation,  h(x,y),  is  known  only  on  a 
discrete  set  of  points,  the  behavior  of  h(x,y)  near  (xp,yp) 
has  to  be  inferred  from  the  values  at  (x  ,y  )  and  at  points 
nearby . 


,»  r.4  ,.l  {.A 


A  planar  approximation  to  the  terrain  elevation  is 
used  in  the  neighborhood  of  the  origin;  that  is 


h(x,y)  =  ax  +  by  +  c 


(3.2-1) 


with  the  parameters  a,  b,  and  c  chosen  to  make  Eq.  3.2-1  the 
best  (least- squares)  planar  fit  to  the  terrain  elevation  di¬ 
rectly  under  the  computation  point  and  the  four  nearest  grid 
points.  Using  the  notation  of  Fig.  3.2-1,  the  values  of  the 
parameters  are 
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Figure  3.2-1  Grid  Points  Used  to  Approximate 

h(x,y)  near  the  Origin 


Note  that,  if  appropriate,  there  is  no  reason  why  a  higher 
order  surface  could  not  be  used. 


-  hl,0  "  h-l,0 
2(Ax) 


b  =  h0,l  ~  h0 , -1 


2(Ay) 


(3.2-2) 


(3.2-3) 


and 


c  = 


h0,0  *  hi>0  1  ho,i  +  h-i,Q  +  ho,-i 

5 


(3.2-4) 


?! 


-a 

k 
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MACSYMA*  (Ref.  7)  was  used  to  obtain  the  integral  of 
Eq.  3.1-2  after  substituting  the  planar  approximation  for  h(x,y) 
in  the  integrand.  The  integration  was  done  by  expanding  the 
integrand  as  a  Taylor  series  in  the  planar  slope ,  and  then  inte¬ 
grating  the  series  term  by  term.  Ideally,  one  would  perform  the 
integration  over  a  rectangle  of  area  A  =  (2Ax)*(2Ay)  centered  at 

the  computation  point.  This  is  not  tractable  in  practice,  how- 

2 

ever,  so  a  circle  of  radius  R  ,  chosen  so  that  nR  =  A,  is  used 

o  o 

instead.  (For  the  DTED  used  in  this  study,  Rq  =  95  m.)  Some 
numerical  experiments  were  conducted  to  determine  how  many  terms 
in  the  series  were  needed  to  obtain  an  accuracy  of  0.01  mgal . 
The  results  showed  that  terms  up  to  degree  2  are  sufficient, 
given  the  ranges  of  terrain  elevations  and  slopes  in  the  DTED 
used  for  this  study.  The  formula  used  to  compute  the  contribu¬ 
tion  to  the  terrain  correction  from  the  area  directly  under  the 
computation  point  is 


IQ(p)  = 


2nGp  <  (R^  + 


2A 

e  )’ 


R. 


e  + 


(■ 


+  6e2R^  +  4e4^  2( 

4(R2  +  e2)*  )  S  i 

°  (3.2-5) 


*MACSYMA  is  a  trademark  of  Symbolics  Inc. 
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2  2  is 

where  s  =  (a  +  b  )’  is  the  slope  of  the  planar  approximation 


and  e  = 


is  the  vertical  distance  from  the  point  P  to 


the  planar  approximation. 

3.3  EXTENT  OF  THE  INTEGRATION  DOMAIN 

The  integration  in  Eq.  3.1-2  cannot  be  performed  numer¬ 
ically  over  the  entire  x-y  plane.  The  domain  of  integration  has 
to  be  truncated,  thereby  introducing  errors  in  the  computation. 
Ideally,  the  truncation  errors  should  be  at  least  as  small  as 
those  caused  by  discretization  effects.  Thus  it  is  appropriate 
to  relate  the  size  of  the  integration  area  explicitly  to  the 
desired  accuracy  of  the  numerical  integration.  Insight  into  the 
dimensions  of  the  required  region  of  integration  is  gained  by 
noting  that  short-wavelength  components  of  gravity  are  due  to 
nearby  masses.  Therefore,  the  integration  in  Eq.  3.1-2  needs 
to  extend  only  far  enough  to  account  accurately  for  these  wave¬ 
lengths  . 

The  error  made  by  limiting  the  integration  in  Eq.  3.1-2 
to  a  circle  of  radius  R  centered  at  the  computation  point  is 


6I(xp,yp,Zp)  =  -Gp 


d 

qr ( r+q) 


dxdy 


(3.3-1) 


In  the  geographical  area  of  interest,  the  values  of  d  are  of  the 
order  of  a  few  hundreds  of  meters.  Since  R  is  expected  to  be 
of  the  order  of  ten  or  twenty  kilometers,  the  approximation 


r  ~  q  is  quite  accurate  for  q>R,  and  Eq.  3.3-1  simplifies  to 


6I(x  ,y  , 
P  P 


V 


=  lG£ 
2 


If 


q>R 


d2  d  d 
—j  dxdy 


(3.3-2) 


Furthermore,  because  it  is  always  true  that  r  >_  q,  the  integrand 
of  Eq.  3.3-1  is  less  than  or  equal  to  the  integrand  of  Eq.  3.3-2; 
therefore,  Eq.  3.3-2  provides  an  upper  bound  on  61. 


Note  that  the  integral  in  Eq .  3.3-2  is  a  convolution, 
which  is  conveniently  rewritten  as 


6i<wv  =  // 


d  yZp ,x,y) 


2  3/2 
f (xp-x)  +  (yp-y)  1 


dxdy 


(3.3-3) 


where  the  integration  domain  is 

A  =  {  ( x , y ) | [ (xp-x) 2  +  (yp-y)2]%  >  R  }  (3.3-4) 


Hence,  the  truncation  errors  can  be  obtained  as  the  convolution 

2 

of  the  squared  height  difference,  d  (z  ,x,y),  and  the  function 


T(x,y) 


if  q(xp,yp,x,y)<R 
if  q(xp,yp,x,y)>R 

(3.3-5) 


If  the  terrain  elevation  is  modeled  as  a  stationary 
random  process,  then  the  squared  elevation  is  also  stationary, 
and  it  follows  from  Eq.  3.3-3  that  the  truncation-error  power 
spectrum  is 


3-6 


(3.3-6) 


*61  =  m2*s 


where  <f>  is  the  power  spectrum  of  the  squared  elevation  dif- 

s  2 

ference,  s(x,y)  =  d  (x,y),  and  T  is  the  Fourier  transform  of 
T  defined  by 


T^i^)  =  JJ  T(x,y)  e 
-00 


-i2n(fjX+f2y) 


dxdy 


(3.3-7) 


This  Fourier  transform  can  be  expressed  as 


T(fi,f2)  =  Q(2nRf ) 


(3.3-8) 


2  2  ie. 

with  f  =  ( f ^ '+  f2)'5.  The  function  Q  is  computed  numerically 


in  terms  of  Jq  ,  the  Bessel  function  of  the  first  kind,  order 


zero 


00 


Q(x) 


f  Jo(u: 

=  J  u2 


du 


(3.3-9) 


For  any  given  choice  of  R,  Eq.  3.3-6  provides  a  decom¬ 
position  of  the  truncation  errors  into  their  spectral  components. 


The  variance  of  the  errors  at  wavelengths  shorter  than  10  km, 
2 

Oq ,  is  obtained  by  integrating  the  spectral  density  over 

radial  frequencies  higher  than  f  =  1/10  cyc/km;  that  is, 


•HJ 


|T|Z4>S  dfxdf2 


(3.3-10) 


f>f 
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Note  that  the  error  variance  oq  depends  on  the  choice  of  R 

through  the  function  T  (see  Eq.  3.3-8).  To  determine  how  large  R 

2 

should  be,  a  table  of  values  of  a  as  a  function  of  R  is  con- 

o 

structed.  The  integration  radius,  R,  is  then  chosen  as  that 
value  in  the  table  for  which  the  error  variance  is  within  the 
required  limits.  Before  such  a  table  can  be  constructed, 
however,  the  spectrum  <|>g  in  Eq.  3.3-10  must  be  evaluated. 

In  the  current  study,  the  spectral  density  of  the 
squared  elevation  differences,  4>g,  was  estimated  using  the  DTED 
in  the  study  area.  Various  east-west  and  north-south  tracks 
were  selected  in  the  area,  and  the  values  of  s  =  d  were 
computed  from  the  terrain  elevation  along  each  track.  One¬ 
dimensional  power  spectra  were  then  obtained  by  fitting  autore- 

2 

gressive  models  to  d  -m  (where  m  is  the  sample  mean  along  the 
track)  using  the  "covariance  method”  (Refs.  8-10).  The  esti¬ 
mated  along- track  spectra  tended  to  cluster  into  two  groups: 
one  for  the  northern,  and  another  for  the  southern  region  of 
the  study  area. 

This  clustering  reflects  differences  in  ruggedness 

which  are  visually  apparent  from  a  graphical  display  of  the 

terrain  elevation  (Fig.  3.3-1).  In  the  range  of  frequencies 
« 

of  interest,  the  spectrum  for  the  northern  tracks  is  uniformly 
higher  than  that  for  the  southern  tracks.  This  behavior  corre¬ 
sponds  to  the  more  rugged  terrain  in  the  northern  region  of  the 
area  being  analyzed. 

It  can  be  seen  from  Eq.  3.3-10  that  the  larger  $s,  the 
larger  the  error  variance  of  the  truncation  errors,  and  the 
more  stringent  the  requirement  on  the  value  of  R.  Therefore, 
the  model  for  the  northern  region  was  used  in  the  determination 
of  the  value  of  R. 


The  function  4>  needed  in  Eq.  3.3-10  is  the  two- 
dimensional  power  spectrum  of  the  squared  terrain  differences. 
Modeling  $s  as  isotropic,  it  can  then  be  computed  from  the  cor¬ 
responding  one-dimensional  along- track  spectrum  using  the  Abel 
transform  (Refs.  11  and  12).  The  result  is  a  spectrum,  <t>(f), 
that  is  a  function  of  radial  frequency  only: 


♦S(f1»f2)  =  *  <f?+f2)1/2 


(3.3-11) 


Since  the  function  T  is  also  only  a  function  of  f  (see 
Eq.  3.3-8),  Eq.  3.3-10  simplifies  to 


|T(f)r<t>(f)  f  df 


(3.3-12) 


Equation  3.3-12  was  used  to  compute  the  data  in  Table  0  3-1, 
which  lists  different  values  of  R  and  the  corresponding  rras 
truncation  error  or  .  For  example,  if  the  vertical  terrain  ef¬ 
fect  is  computed  for  wavelenghts  shorter  than  10  km  using  DTED 
within  a  3-km  radius  of  the  computation  point,  then  Table  3.3-1 
indicates  that  the  rms  error  due  to  the  finite  amount  of  DTED 
is  0.0081  mgal . 
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4. 


DISCUSSION 


4 . 1  DATA 

The  Defense  Mapping  Agency  provided  the  DTED  used  for 
this  analysis.  These  DTED  cover  a  region  that  includes  the 
Clinton-Sherman  Airport  in  Oklahoma  and  consist  of  terrain 
heights  above  sea  level  on  a  uniform  latitude-longitude  grid 
having  a  point  spacing  of  three  seconds  of  arc,  with  the  heights 
expressed  to  the  nearest  meter.  The  object  of  the  present 
analysis  is  to  estimate  the  high-frequency  vertical  gravity 
field  in  a  512  by  512-point  subregion  that  includes  a  topo¬ 
graphic  feature  called  Baker  Peak.  It  is  desired  to  compute 
the  local  terrain  effect  at  each  grid  point  in  this  subregion 
with  an  rms  accuracy  of  0.01  mgal  for  wavelengths  shorter  than 
10  km.  According  to  Table  3.3-1,  to  achieve  this  accuracy, 
DTED  are  needed  within  a  3-km  radius  of  each  computation  point. 

Therefore,  the  numerical  quadrature  of  the  terrain-correction 

2 

integral  was  performed  over  a  square  area  spanning  9 n  km  and 
centered  on  each  computation  point. 


In  Fig.  3.3-1,  the  DTED  are  rendered  as  a  three-dimensional 
surface,  with  the  vertical  scale  greatly  magnified  for  visual 
clarity.  Although  the  northern  (upper)  half  of  Fig.  3.3-1  looks 
mountainous,  the  DTED  actually  fall  in  a  range  from  approxi¬ 
mately  200  m  to  700  m  above  sea  level.  The  region  depicted  in 
Fig.  3.3-1  extends  57  km  east-west  and  65  km  north-south.  The 
average  distance  between  grid  points  is  76.5  m  east-west 
and  92.7  m  north- south. 
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4.2 


COMPUTATIONAL  APPROACH 


To  generate  estimates  of  high-frequency  gravity  on  a 
dense  grid,  the  local  terrain  effect  is  first  computed  at  each 
grid  point.  Then  the  resulting  two-dimensional  (2-D)  field 
is  filtered  with  a  2-D  high-pass  filter.  Computing  the  terrain 
effect  is  the  major  data  processing  effort,  because  the  terrain 
effect  integral  in  Eq.  3.1-2  must  be  evaluated  at  each  point 
on  the  grid.  In  particular,  for  the  512  by  512-point  grid  used 
in  this  paper,  the  terrain  effect  integral  is  evaluated 
numerically  262,144  times.  Each  of  these  numerical  quadratures 
involves  the  following  three  steps:  (1)  the  integrand  is 

evaluated  at  the  4,189  grid  points  within  the  54  by  71  point 
domain  of  integration;  (2)  at  each  of  these  points,  the  value 
of  the  integrand  is  multiplied  by  a  weighting  factor,  w^, 
determined  by  the  numerical  quadrature  rule;  and  (3)  the  re¬ 
sulting  products  are  summed  and  scaled  by  the  grid- spacing 
factor  4-Ax*Ay. 

The  computational  approach  described  here  is  appropri¬ 
ate  for  computers  having  a  Cray-type  architecture  for  parallel 
vector  computation,  e.g.,  the  Alliant  FX/8  minisupercomputer 
used  in  this  study.  On  such  machines,  arithmetic  operations 
applied  uniformly  to  large  arrays  can  be  performed  much  more 
rapidly  than  an  equivalent  number  of  scalar  calculations. 
Therefore,  to  achieve  high  computational  speed,  the  computations 
are  expressed  as  uniform  operations  on  large  arrays. 

The  terrain  effect  integral  is  evaluated  using  the 
following  nine-point  numerical  quadrature  formula  (Ref.  13): 


//  f (x ,y )  dxdy  =  4*Ax*Ay  ^  f^xk,yk^wk  +  (4.2-1) 


where  the  domain  S  is  the  2*Ax  by  2* Ay  rectangle  (-Ax  <  x  <_  Ax, 
-Ay  <  y  <  Ay)  of  area  A  =  4*Ax*Ay.  As  depicted  in  Fig.  4.2-1, 
the  knot  points  (x^y^)  are  on  a  rectangular  grid  with  spacings 
Ax  and  Ay  in  the  x  and  y  directions,  respectively,  and  the 
weights,  4w^,  are  the  outer  product  of  the  weights  (1/6,  2/3, 
1/6)  for  the  three- term  Simpson's  formula.  The  error  of  this 


quadrature  rule  is  asymptotically  of  order  A  as  A  goes  to 


zero. 


According  to  the  error  analysis  in  Section  3.3,  the 
domain  of  the  terrain-effect  integral  can  be  truncated  to 
a  rectangular  area  centered  on  the  computation  point  P  as  de¬ 
picted  schematically  in  Fig.  4.2-2.  This  domain  contains 
(2M+1) * (2N+1)  DTED  grid  points.  Applying  Eq.  4.2-1  to  each 
adjacent  9-point  rectangle  and  summing  the  results  over  the 
domain  depicted  in  Fig.  4.2-2  will  yield  the  desired  numerical 
quadrature,  except  for  the  center  region  containing  point  P, 
which  is  omitted  from  this  quadrature.  The  center  region  is 
handled  separately  as  discussed  in  Section  3.2.  To  achieve 
high  speed  on  a  parallel-vector  computer,  the  multiple  appli¬ 
cations  of  Eq.  4.2-1  should  be  formulated  as  a  single  large 
array  operation.  This  is  done  by  developing  a  single  (2M+1) 
by  (2N+l)-point  (large-scale)  quardrature  formula  that  is  equiv¬ 
alent  to  multiple  applications  of  the  9-point  formula.  The 
resulting  weight  pattern  over  the  interior  of  the  domain  is 
depicted  in  Fig.  4.2-2.  (The  weights  along  the  outer  edges 
and  the  weights  along  the  omitted  inner  region  containing  point 
P  are  different  from  the  interior  pattern  depicted  in  Fig.  4.2-1) 


After  evaluating  the  terrain  correction  integral  on 
the  512  by  512  grid,  the  final  step  in  estimating  high-frequency 
gravity  is  to  filter-out  the  low-frequency  components  of  the 
computed  terrain  effects.  A  2-D  finite-impulse-response  high- 
pass  filter  is  used  for  this  purpose.  The  output  of  the  filter, 
y( j ,k) ,  is  the  value  of  the  input  datum,  x(j,k),  minus  the 
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Figure  4.2-1  Domain  and  Weights  W^  of  Nine-Point 

Quadrature  Formula 
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Schematic  Representation  of  Integration 
Domain  for  Point  and  Interior  Weight  Pattern 
for  Large-Scale  Numerical  Quadrature 


average  of  x(j,k)  over  a  (2N  +1)  by  (2N  +l)-point  rectangl* 

x  y 

centered  on  point  (j,k): 


Ny  Nx 

y(j,k)  H(m>n)x(j-m>k-n) 

m=-Ny  n=-Nx 

H(0 ,0)  =  1  -  ( 2N  +1 ) " 1 ( 2N  +1 ) " 1 

x  y 

H(m,n)  =  - ( 2N  +1 ) " 1 ( 2N  +1 ) " 1 
x  y 


(4.2-2) 


(4.2-3) 


(4.2-4) 


The  transfer  function  of  this  filter,  H(Fy,Fx),  is  expressed 

using  the  normalized  frequencies,  Fx  =  Axfx  and  Fy  =  Ayfy, 

where  F  and  F  range  in  absolute  value  from  0  to  1/2: 
x  y 


H(Fy,Fx)  =  1  -  Hy(Fy)Hx(Fx) 


Hx(Fx}  =  sinc[(2Nx+l)Fx]/sinc(Fx) 


Hy(Fy)  =  sinc[ (2Ny+l)Fy]/sinc(Fy) 


sinc(0)  =  1 


sinc(x)  =  sin(nx)/(nx) 


for  x  /  0 


(4.2-5) 


(4.2-6) 


(4.2-7) 


(4.2-8) 


(4.2-9) 


The  frequency  response  of  this  filter  (i.e.  ,  the  squared  magni¬ 
tude  of  the  transfer  function)  is  plotted  in  Fig.  4.2-3  for 
F  =  0  as  a  function  of  the  dimensionless  east  frequency  w  f  , 

y  XX 

where  wx  is  the  length  of  the  filter  impulse  response  in  the  x 
direction,  and  fx  is  the  frequency  expressed  in  cycles  per  unit 
length  (i.e.,  wxfx  =  (2Nx+l)Axfx  =  (2Nx+l)Fx).  The  low-frequency 
limit  of  this  filter  is  defined  as  the  frequency  for  which  the 
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Figure  4.2-3  Frequency  Response  of  High-Pass  Filter 


for  East  Frequencies  ( f ^  =  0) 


squared  mangitude  of  the  transfer  function  is  1/2.  It  can  be 


shown  that  a  specified  low-frequency  limit,  f  q  ,  is  achieved  for 


both  north  and  east  frequencies  by  selecting  the  filter  para¬ 


meters  N  and  N  to  satisfy  the  following  formulas,  in  which 
x  y 


int(x)  denotes  the  value  of  x  rounded  to  the  nearest  integer: 


_  .  JO. 76  l\ 

-  int(2Axf0-  2) 


(4.2-10) 


.  .  JO. 76  l\ 

-  int(2Ayf0-  2) 


(4.2-11) 
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Figure  4.3-1 


Three-Dimensional  Rendering  of 
Estimated  Vertical  Gravity  for 
Wavelengths  Shorter  than  10  km 


4.3  RESULTS 


The  results  of  using  the  computational  approach  des¬ 
cribed  in  Section  4.2  are  depicted  in  Fig.  4.3-1,  which  is  a 
three-dimensional  rendering  of  the  estimated  vertical  gravity 
field  for  wavelengths  shorter  than  10  km.  The  elevation,  zp , 
is  equal  to  the  average  of  the  DTED  over  the  analyzed  area. 
The  northern  (upper)  region  in  Fig.  4.3-1,  which  contains 
1.2  x  105  grid  points,  has  an  rms  value  that  is  five  times  the 


rms  field  in  the  southern  region,  which  contains  1.2  x  10  grid 
points.  Key  statistics  for  the  estimated  gravity  fields  in 
these  two  regions  are  compared  in  Table  4.3-1. 


The  error  analysis  in  Section  3.3  indicates  that  the 
high-frequency  vertical  terrain  effect  was  computed  with  an 
rms  error  of  0.0081  mgal,  because  the  DTED  data  extended  3  km 
from  each  computation  point.  To  verify  this  small  error,  the 
terrain  effects  were  recomputed  using  DTED  data  within  a  10-km 
radius  of  each  computation  point.  According  to  Table  4.3-1, 
the  rms  error  of  this  computation  is  0.00054  mgal.  This  more 
accurate  field  calculation  was  subtracted  from  the  original  re¬ 
sults,  and  the  rms  difference  was  consistent  with  the  theoreti¬ 
cal  rms  error  of  0.0081  mgal  given  in  Table  3.3-1  for  R  =  3  km. 


An  additional  error  analysis,  not  presented  here, 
shows  that  a  much  larger  extent  of  DTED  (i.e.,  larger  R) 
is  required  to  achieve  similar  accuracies  when  computing  the 
deflections  of  the  vertical  contributed  by  the  terrain  effect. 
For  example,  an  rms  accuracy  of  0.15  mgal  requires  DTED  within 
a  10-km  radius  of  each  computation  point. 


TABLE  4.3-1 


STATISTICS  OF  ESTIMATED  GRAVITY  FIELDS 
FOR  WAVELENGTHS  SHORTER  THAN  10  km 
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V 

v’ 

GRAVITY  IN 

GRAVITY  IN 

*«' 

M 

STATISTIC 

NORTHERN  REGION 

SOUTHERN  REGION 

(mgal) 

(mgal) 

1 

MAXIMUM 

18.0 

1.5 

ft 

MINIMUM 

-6.8 

-2.2 

STD .  DEV . 

3.0  c 

0.6  , 

NO.  POINTS 

1.2x10° 

1.2xl04 

1 

4.4 


SUMMARY  AND  APPROPRIATE  NEXT  STEP 


The  formulation  developed  in  this  report  surmounts  cer¬ 
tain  mathematical  difficulties  formerly  associated  with  cal¬ 
culating  that  portion  of  the  gravity  field  which  is  induced  by 
visible  terrain.  It  also  uses  to  advantage  modern  parallel¬ 
processing  computer  architectures  to  optimize  calculations 
which  would  be  arduous  or  impractical  if  performed  in  scalar 
sequential  fashion.  A  statistical  error  analysis  provides  the 
means  to  control  accuracy.  The  approach  is  applied  to  digital 
terrain  elevation  data  from  the  Clinton-Sherman  test  range  to 
develop  both  a  set  of  terrain-effect  gravity  values  and  an  under¬ 
standing  of  the  accuracy  with  which  the  field  can  be  determined. 

The  accuracy  uncertainty  associated  with  the  DTED 
motivates  further  efforts  to  understand  and  quantify  the  data 
errors.  Is  there  a  frequency  range  over  which  the  accuracy 
does  approach  its  precision?  Can  the  errors  be  predicted? 

Can  they  be  modeled?  Resolution  of  these  issues  could  provide 
the  impetus  for  turning  DMA's  DTED  base  into  a  very  useful 
tool  for  supplementing  gravity  data  needed  to  support  DoD 
weapon  systems  programs. 
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